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ANALYSIS AND A NINE NODE LINEAR CRACK ELEMENT 
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I SUMMARY 

This paper shows how a two-dimensional crack element was implemented into 
NASTRAN as a user dummy element and used to study failsafe characteristics of the 
[C5A fuselage. The element is formulated from Reissner's functional requiring that it 
satisfy compatability with the linear boundary displacement elements in NASTRAN. 
Its accuracy is demonstrated by analyzing for the stress intensity factors of two simple 
[crack configurations for which there are classic solutions. 


1 INTRODUCTION 

A Lockheed requirement during the design of the C5A was for the pressurized 
fuselage structure to have the capability of sustaining a longitudinal crack in the skin. 
Circumferential straps are attached to the skin midway between the frames to provide 
this capability. Conventional methods for designing the straps and the system for 
attaching them to the skin are conservative and as a result dictate the use of high strength 
fasteners. It is now possible, due to the introduction of special crack tip singularity 
elements, to make a finite element analysis which eliminates much of this conservatism. 

I 

Part of the criterion for determining if a crack can be arrested is to find a 
f configuration for which the stress intensity factor at the tip of the crack is lower than 
the critical stress intensity factor for the skin material. NASTRAN was considered the 
best analysis tool with which to analyze the various crack configurations. This was due 
to the ease with which a crack element can be incorporated into it as a user dummy 
element. 
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LIST OF SYMBOLS 


Stress 

Strain 

Displacement 

Vector normal to a surface 
Applied tractions 
Applied displacements 
Complementary energy function 
Body forces 
Volume of an element 

tj 

Surface on which T* are specified 

Surface on which U* are specified 

Interior boundary subdividing a region into finite elements 

Compliance matrix 

Shear modulus 


Poisson's ratio 


SINGULARITY CRACK ELEMENT FORMULATION 


In this section the stiffness matrix and the relationship between the nodal 
displacements and the stress intensity factors for a crack element are developed. 

Reissner's functional. Reference 1, fora region subdivided by boundaries D into 
finite elements gives 



where the summation includes all the finite elements making up the region 

If the displacements of the surfaces D are assumed compatible between finite 
elements, then the necessary conditions for J R to be stationary due to a variation 


in the functions a lK , and 

are the following five Euler equations: 
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A modified form of (1) is obtained by neglecting body forces and restricting 
the functions a LK and to the set which satisfy (3), (4), and (5) 


( 8 ) 


'.-E 


-/w c (<7 l ,)dv+ j, 

Jv Js t 




°lk VtcVfdS + a lK rj K H L dS 


For this form of the functional the displacements // t only appear in the boundary 
integrals and hence do not have to be defined in the interior of the element. 


In applying J, to the finite element shown in Figure 1 the stresses a LK are 
approximated by finite series and the displacements fi L are required to vary linearly 
between nodal displacements { Q } . This ensures displacement compatibility with 
the linear boundary displacement elements already in the NASTRAN system. Expressing 
these approximate functions in matrix form 


(9) 

\° ! = 

[p] 1*1 

stresses in V 

(10) 

\t*o} - 

L l »] \°\ 

displacements of the surface D 

(ID 

f T 0 } = 

[p»] {*1 

o LK rj K at the surface D 


where £ P 1 are the assumed series of stress functions and are the corresponding 

coefficients, are linear interpolations of the nodal displacements j Q j along 

the boundary D, and £ P D J are the values of |_ P ]] that give the stresses | a j in 
Equation 9 normal to the boundary D. The stresses on the boundary Sj, do not have 
to be defined as there are no enforced boundary displacements associated with this 
element. 
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to J . 


Substituting (9), (10), and (11) into (8) gives the contribution of the crack element 


(12) j, c = - [kc]| 0 | + |_/»J[ka]{q ( 


where 


(13) [KC] = / [ P ] T [ S ][ P ] dV 

Jy 

(14) [KA]= f [n][*Jis 


Due to the fact that the coefficients , { 0 } , are independent between elements the 
necessary condition for J lc to be stationary is dJ] C q . This yields: 

w x 

(15) & = [kc] [ka] { Q ( 


Substituting (15) into (12) and finding a stationary value of Jjq with respect to | Q | 
give the contribution of the crack element to the total stiffness matrix 



Each term in the series used to approximate the stresses has to satisfy equilibrium 
and the stress boundary conditions. The series also has to include the known stress 
singularity near the tip of the crack. The form that has been assumed for this element 
is that which was developed from the Williams series of stress functions by Aberson and 
Anderson in Reference 2. The two dimensional stress distributions in polar coordinates 
are: 
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(17) a„(R,0) ='^0sJ n / 4 R /2 1 [-(n + 2)Cos(n/ 2+ l)0 + f( n )(n-6)Cos( n / 2 -l)0 


2>p 


% R /2 1 [g(n)(n + l)Sin( n /2+l)0-(n-6)Sin(n/ 2 -l)0 


(18) oAR,9) 
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%-ir 

'/ 4 R (n + 2)Cos( n / 2 + l)d-f(n)Cos( n /2-l)0(n+2) 




n/,-1 


!4 R r9(n)(n + 2)Sin(% + 1 )# + ( n + 2)Sin(%-l )0 


T r ^( R , 0 ) = y^ £S n jn / 4 R /2 1 (n + 2)Sin( n / 2 + l)0-f(n)(n-2)Sin( n / 2 -l)0 
' n=i ^ L 


•E " An " 


n/, - 1 

'/ 4 R g(n)(n+2)Cos("/ 2 +l)0-(n-2)Cos("/ 2 -l)0 


in which 


n /i + 1 

n / 2 +H) n 


g(n) = 


n / 2 -(-ir 
n /2+ 1 


where /3S n and /SA n are the symmetric and antisymmetric parts of j 
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The symmetric and antisymmetric coefficients and 0A, are associated with 
the singularity R /2 . They are related to the opening and sliding mode stress 

intensity factors K, and K n by the following formulas 


(20) K, = 3 V2/7 0S, 


1(21) K„ = Vrfl 


Hence once the displacements j Q } have been determined the leading coefficients 
/^cind 0S, can be recovered through (15) and subsequently the stress intensity 
factors through (20) and (21). 


! IMPLEMENTING THE CRACK ELEMENT INTO NASTRAN 

The nine node crack element has been implemented into NASTRAN through the 
dummy element capability. This capability permits a user to enter his own element 
subroutines for the purpose of generating the stiffness and mass matrix contributions, 
the thermal load contributions and for the computation of various stresses and forces 
for output (see section 8.8.5 of the NASTRAN PROGRAMMERS MANUAL, Reference 
3). This procedure is relatively simple compared to adding an entirely new element 
to the system. 

The crack element has been implemented as a DUM2 element. The format for the 
: ADUM2, CDUM2, and PDUM2 bulk data cards which are used to enter the geometry, 
property, and connectivity data is shown in Figure 2. The procedure used to 
implement the element is as follows. 

o Create an element stiffness subroutine KDUM2 which computes and outputs 
to functional module SMA1 one 6x6 matrix for each connecting grid point 
with respect to the connective pivot point, 

o Create two subroutines SDUM21 and SDUM22 to compute and output to 
functional module SDR2 the stress intensity factors. 
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o 


Remap LINK3 to include the new routine KDUM2 and LINK13 to include 
the new routines SDUM21 and SDUM22. 


All the element stiffness subroutines called by SMA1, including KDUM2, are 
overlayed. Therefore to avoid reducing the working core available to SMA1, KDUM2 
was programmed in less core than the largest amount used by any of the existing element 
stiffness routines. To do this it was necessary to fix the shape of the element shown in 
Figure 1. Chosen was a square with the nine grid points equally spaced around the 
boundary. This made it possible to compute the integrations involved in (13) 

and (14) externally to NASTRAN for a unit element size. Further (13) was 
integrated for a unit value of each of the two independent coefficients for an isotropic 
material in the compliance matrix fs 1 i.e. 
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S, 
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1 / 2 G 

-U/2G PLANE STRAIN 

-u/2G(l + l>) PLANE STRESS 


The resulting unit matrix for [KAj and the two unit matrices for 
permanent data in the subroutines KDUM2 and SDUM21. 



are used as 


Subroutine KDUM2 forms the matrices [ka] and [kcJ in double precision from 
the three unit matrices for a specific element size and material. It then uses these in 
(16) to compute the 2x2 stiffness matrices in the local element coordinate 
system for each grid point associated with the respective connective pivot point. Finally 
it transforms these 2x2 matrices into 6x6 matrices in the global coordinate system. 


188 


i 


Similarly subroutine SDUM21 forms the matrices [kaJ and in single 

precision for a specific element size and material. It then computes the two stress 
intensity factors K, and for a unit value of each of the grid point displacements 

{ Q } in the global system using 15, 20, and 21. Subroutine SDUM22 
| computes the final stress intensity factors for specific grid point displacements. 

i 

! 


| ELEMENT EVALUATION 

The capability of the element to predict accurate stress intensity factors has been 
demonstrated by numerous analyses. Two are presented here for which there are 
| known classic solutions. 

The opening mode stress intensity factor for an isotropic materia! can be expressed 
in the form. Reference 4. 


I (22) K 1 = oVna f 


Where a \ZJfa is the stress intensity factor fora central crack of length 2a 
in an infinite plate loaded by a far field stress a acting normal to the crack, f 
is either the correction factor associated with Bowie's analysis for the presence of a 
hole. Reference 5, or Isida's analysis for a finite width plate. Reference 6. 

I 

The finite element model shown in Figure 3 represents a crack in a finite width 
plate loaded by a far field stress acting normal to the crack. The plate is 81 cm 
wide 168 cm long and the crack length is varied between 10 and 51 cm. The 
model which idealizes one quarter of the plate represents the full structure thru the use 
of symmetric boundary constraints. The crack element which overlaps a symmetry 
boundary is forced to deflect symmetrically thru the use of multi-point constraint 
equations. It should be noted that when the multipoint constraint equations are used in 
this way the element that is being forced to act symmetrically should be specified with 
half its actual thickness. Besides the DUM2 crack element the model consists of 66 
CQDMEM elements, and 36 CTRMEM elements. It has 124 grid points, 219 active 
freedoms, and a semiband width of 37. Figure 3 shows that for this type of 
idealization the crack element computes the stress intensity K, to within 3% of that 
given by equation (22). 
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The finite element model shown in Figure 4 represents symmetric cracks protruding 
from a hole in the center of an infinite plate. The plate is 102 cm wide and 102 cm 
long. The hole has a diameter of 5 cm and the crack lengths are varied 
between 0.8 and 8 cm . The model uses the same idealization techniques employed 
in the previous example to represent the plate. It consists of 154 CQDMEM elements, 
5 CTRMEM elements and one DUM2 crack element. Figure 4 shows the results to be 
within 3% of equation 22. 


C5A FAILSAFE ANALYSIS 

The C5A fuselage has a failsafe criterion which requires that a 30 cm longitudinal 
crack in the cover skin will not result in a catastrophic failure when the structure is 
subjected to a normal operating internal pressure. To satisfy this requirement, 
circumferential failsafe straps are attached to the skin mid way between the frames for 
the purpose of arresting such a crack, see Figure 5. An analysis based on Lockheed 
data sheets is conservative and as a result dictates the use of expensive high strength 
fasteners to attach the strap to the skin. The following analysis using NASTRAN shows 
that less expensive aluminum rivets can be used instead. 

As a longitudinal skin crack passes under the failsafe strap the stress intensity at 
the tip reduces. The crack will cease to propagate if the stress intensity becomes less 
than the critical stress intensity for the material, provided that the fasteners in the 
strap do not fail first. The finite element model shown in Figure 6 represents a typical 
region of the C5A aft fuselage. It considers the frame at fuselage station (F.S.) 1804 
to be failed and is used to analyze various lengths of a skin crack which propagates 
towards the failsafe straps at F.S. 's 1794 and 1814. The model is two-dimensional , i .e. 
fuselage curvature and out of plane deflections are ignored. Advantage has been taken 
of two symmetry planes by idealizing only one quarter of the actual damaged region. 

The crack element which lies across a symmetry boundary is again forced to displace 
symmetrically through the use of multi-point constraint equations. The frame cap and 
the skin are represented as an integral structure through the use of CROD,CTRMEM, 
CQUAD, and the CDUM2 elements. The strap which is considered as a separate unit 
uses the CTRMEM and CQUAD elements. The twelve fasteners closest to the crack are 
each idealized by a system of CROD elements. The remaining fasteners are lumped into 
groups of approximately 4 and idealized by CELASI elements. The model is loaded by 
concentrated forces which represent the hoop loading on a 244 cm radius structure for a 
normal operating internal pressure with a dynamic factor of 1 . 15. The loading is 
applied to the model in proportion to the circumferential cross sectional area. This 
does not represent the true circumferential loading as it neglects the bulging effect of 
the skin between the frames. It is conservative however in that it overloads the frame 
at the center of the crack. 
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i In order to show the fasteners capable of carrying the transfer load it is 
^ecessary to consider the nonlinear load deflection response for the twelve fas- 
'teners closest to the crack, see figure 7. Because the CELASI elements do not 
have the capability of representing nonlinearity it was necessary to idealize 
'these fasteners by the system of CROD elements shown in Figure 8. The rods con- 
nect the coincidental grid points A and B on the skin and strap respectively 
through the grid point C. The elements have a combined length of 2.54 cm (1 in.) 
and a cross-sectional area of 6.45 cm2 (i i n 2). The use of English units allows 
the load deflection curves in Figure 7 to be input directly on Tables 1 cards as 
a stress strain curve for the rod elements. The only other region of the struc- 
ture to experience plasticity is the tip of the crack and since this is ignored 
in linear fracture mechanics the most economical way of executing the analysis 
is to divide the model into two substructures: the skin, frame, and strap being 

j included in substructure 1 and the fasteners in substructure 2. Substructure 1 
is analyzed first for an increment of the external loads using rigid format 1. 
jThe freedoms for the grid points common to the fasteners are included in the 
'A' set and the rigid format altered to terminate once the reduced A set stiff- 
ness and loads matrices are formed. Substructure 2 is then analyzed using rigid 
[format 6, the piecewise linear analysis. Rigid format 6 is altered to read the 
I A set stiffness and loads matrices for substructure 1 and add them to the appro- 
priate terms in the G set stiffness and loads matrices for substructure 2. The 
results of the piecewise linear analysis give the desired fastener loads. To 
obtain the stress intensity factor it is necessary to restart the analysis for 
! substructure 1 using the A set displacements resulting from the analysis of 
substructure 2. 

Substructure 1 consists of one CDUM2, 56 CROD, 71 CTRMEM and 814 CQDMEM 
elements. It has 952 grid points, 1876 active freedoms and a semiband width 
of 84. There are 136 freedoms in the A set for which it took 31 CPU minutes, 
on a UNIVAC 1106 computer executed in a time sharing mode, to form the reduced 
stiffness and loads matrices. Substructure 2 consists of 56 CELASI elements and 
24 CROD elements. It has 80 grid points, 148 active freedoms and a semi band 
width of 142. It took 10 iterations to obtain the elastic-plastic solution 
using 14 CPU minutes. The back feed into substructure 1 to obtain the stress 
intensity factor ran for 8 CPU minutes. All the above run times are for 
Level 15.1. 

The analysis has been made for both the original 3.2 mm (1/8 in.) 

YBO TAPERLOKS and the proposed 4.0 mm (5/32 in.) aluminum rivets. Four crack 
lengths were considered for each system, the results being plotted in Figures 9 
and 10 respectively. They show the ratio of the stress intensity to the critical 
stress intensity (K]/K c ) and the ratio of the maximum fastener load to the allow- 
able fastener load ( P/P a ) plotted against the half crack length. Figure 9 shows 
that when the TAPERLOKS are used the stress intensity reduces to the critical 
value at a half crack length of 22 cm. At this length the fasteners are only 
working to 40% of their allowable load; hence they are shown to provide the 
necessary failsafe characteristics. This is only to be expected as the more 
conservative Lockheed data sheets also show the TAPERLOKS capable of arresting 
the crack. Figure 10 shows that when the rivets are used the stress intensity 
reduces to the critical value at a half crack length of 24 cm. At this length 
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the rivets are working to 67% of their allowable load therefore, unlike the con- 
ventional analysis, this analysis shows that they also are capable of providing 1 
the necessary failsafe characteristics. This analysis was validated by a test 
program which demonstrated that the rivets are capable of arresting the skin 
crack. 


CONCLUSION 

The NASTRAN capability of allowing a user to implement a dummy element was 
found to be relatively simple to use and exceedingly useful. Lockheed Georgia 
Company has plans to implement more elements into the system. In particular work is 
underway on a fastener element, similar to CELAS1, that has the capability of 
representing a nonlinear load deflection curve. 


REFERENCES 

1. Fung, Y.C.: Foundations of Solid Mechanics. Prentice-Hall Inc., 1965, pg 300. 

2. Aberson, J.A., Anderson, J.M.: Cracked Finite Elements Proposed for NASTRAN. 
NASTRAN: User's Experiences. NASA TM X-2893, 1973 

3. The NASTRAN Programmer's Manual . NAS SP-223 (01), Sept. 1972. 

4. Bowie, O. L.: Analysis of an Infinite Plate Containing Radial Cracks Originating 
from the Boundary of an Internal Circular Hole. Journal of Mathematics on 
Physics, Vol. 35, 1956. 

5. Isida, M.: On the Tension of a Strip with a Central Elliptical Hole. Transactions, 
Japanese Society of Mechanical Engineers, Vol. 22, 1956. 

6. Wilhem, D.P.: Fracture Mechanics Guidelines for Aircraft Structural 
Application. AFFDL-TR-69-1 1 1 , February 1970. 


192 




PDUM2I 

PDUM2 


PIP 

1 


MID 


N 

1 


WHERE 

NG 

NC 

NP 

ND 

EID 

PID 

Gi 

MID 

t 

N 


NUMBER OF GRID POINTS CONNECTED BY DUM2 CRACK ELEMENT 

NUMBER OF ADDITIONAL ENTRIES ON CDUM2 CONNECTION CARD 

NUMBER OF ADDITIONAL ENTRIES ON PDUM2 PROPERTY CARD 

NUMBER OF DISPLACEMENT COMPONENTS AT EACH GRID POINT 

ELEMENT IDENTIFICATION NUMBER 

IDENTIFICATION NUMBER OF A PDUM2 PROPERTY CARD 

GRID POINT IDENTIFICATION NUMBERS OF CONNECTION POINTS 

MATERIAL IDENTIFICATION NUMBER 
ELEMENT THICKNESS 
N=1 INDICATES PLANE STRESS 
N=0 INDICATES PLANE STRAIN 


Figure 2. Bulk Data Cards For a DUM2 Crack Element 


I9*t 



































1 



Figure 4. Ki Comparison For Symmetric Cracks at a Hole In an Infinite Plate 
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SYMMETRIC BOUNDARY 





Figure 6« NASTRAN Model of a C5A Fuselage Skin Crack 
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SYMMETRIC BOUNDARY 








Ki/Kc or P/Pa Ki/Kc or P/Pa 



HALF CRACK LENGTH a/2, cm 

Figure 10. NASTRAN Results From the Finite Element Model 
Shown in Figure 6 for 4.0 mm Aluminum Rivets. 




